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ABSTRACT 

In classical loss cone theory, stars are supplied to a central black hole via 
gravitational scattering onto low angular momentum orbits. Higher feeding rates 
are possible if the gravitational potential near the black hole is non-axisymmetric 
and the orbits are chaotic. Motivated by recently published, self-consistent mod- 
els, we evaluate rates of stellar capture and disruption in triaxial nuclei. Rates 
are found to substantially exceed those in collisionally-resupplied loss cones, as 
long as an appreciable fraction of the orbits are centrophilic. The mass captured 
by a black hole after a given time in a steep (p ~ r~ 2 ) nucleus scales as <r 5 with 
a the stellar velocity dispersion, and the accumulated mass in 10 10 yr is of the 
correct order to reproduce the M» — a relation. Triaxiality can solve the "final 
parsec problem" of decaying black hole binaries by increasing the flux of stars 
into the binary's loss cone. 

1. Introduction 

Fueling of active galactic nuclei (AGNs) or quasars requires accretion at rates of M « 
1.5e^jL 46 A / ioyr _1 onto the central supermassive black hole (SBH), where L 46 is the energy 
output in units of 10 46 ergs s _1 and eo.i is the mass conversion efficiency in units of its 
canonical value 0.1. The fueling problem is usually broken into three parts: what is the 
fuel; how is it channeled into the SBH; and how does it radiate a substantial fraction of 
its energy before disappearing down the hole? This paper addresses the first two questions. 
Galactic spheroids have ample supplies of both stars and gas, but it is difficult to come up 
with mechanisms that can extract almost all of a mass element's orbital angular momentum 
in a few crossing times, as required if the matter is to find its way into the event horizon of 
the black hole. Gas release through tidal disruption of stars is a possible fueling mechanism 
(Hills 1975; Frank & Rees 1976), but the rate at which stars are scattered onto low angular 
momentum orbits by gravitational encounters is too low to reproduce observed luminosities 
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( Young, Shields & Wheeler 1977; Frank 1978). The fact that most AGNs are in spiral galaxies 
suggests that interstellar gas provides the bulk of the fuel. Gas can be driven into the SBH 
by torques from non-axisymmetric potential perturbations, due to stellar bars or to transient 
distortions of the potential during mergers or accretion events (Shlosman, Begelman & Frank 
1990). 

As a number of authors have pointed out, stellar feeding rates are also enhanced by the 
presence of non-axisymmetric perturbations (Norman & Silk 1983; Gerhard & Binney 1985). 
A barlike or triaxial potential is populated mainly by box orbits, and a star on a box orbit 
passes near the center once per crossing time. If a SBH is added, most of these "centrophilic" 
orbits become chaotic due to large-angle deflections by the SBH (Valluri & Merritt 1998). 
Estimates of feeding rates due to centrophilic orbits in triaxial potentials are several times 
larger than rates due to scattering onto eccentric orbits in the spherical or axisymmetric 
geometries (Gerhard & Binney 1985). While such an enhancement is significant, it is still 
not enough to explain the high luminosities of AGNs and quasars. 

A number of factors motivated us to re-open the question of stellar fueling rates in non- 
axisymmetric nuclei. The M, — a relation demonstrates a tight link between SBH masses 
and the kinematics of their stellar spheroids. While the origin of the relation is still uncer- 
tain, some scenarios postulate a role for stellar feeding (e. g. Zhao, Haehnelt & Rees 2002). 
Capture of stellar- mass objects by SBHs may provide an important source of signals for grav- 
itational wave detectors like LISA (Hughes et al. 2001). Imaging of the centers of galaxies on 
scales of ~ 10 pc reveals a wealth of features in the stellar distribution that are not consistent 
with axisymmetry, including bars, bars-within-bars, nuclear spirals, and other misaligned 
structures (Wozniak et al. 1995; Rest et al. 2001; Peng et al. 2002; Erwin & Sparke 2002). 
The steep power-law dependence of stellar density on radius near the centers of many galax- 
ies revealed by HST (Crane et al. 1993; Gebhardt et al. 1996) implies higher stellar feeding 
rates than in older models (Norman & Silk 1983; Gerhard & Binney 1985) which postulated 
constant-density cores. Triaxiality has recently been shown to be sustainable in numerical 
models for black-hole nuclei (Poon & Merritt 2001, 2002, 2003, hereafter Papers I-III). Fur- 
thermore such models can contain a large population of centrophilic - typically chaotic - 
orbits, as high as ~ 75%. 

This paper, the fourth in a series on the dynamics of triaxial black-hole nuclei, examines 
the behavior of centrophilic orbits and the implications for black hole feeding. Test-particle 
integrations in fixed potentials are combined with knowledge of the orbital population in 
the self-consistent models to infer the rate at which stars would be supplied to the central 
SBH. We confirm and extend the results of earlier authors, who investigated stellar feeding 
rates in triaxial nuclei with constant-density cores. We find that stellar capture rates in 
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steep power-law nuclei can be orders of magnitude greater than in nuclei with cores; in fact, 
feeding rates can meet or even exceed the so-called "full loss cone" rate, i.e. the capture rate 
in a spherical or axisymmetric galaxy in which the loss cone is continuously repopulated. 
The accretion rate in a dense, p ~ r~ 2 nucleus scales as <r 5 and the mass accumulated in 
10 10 yr is of the right order to reproduce the M, — o relation. 

A number of consequences follow from such high rates of stellar feeding. Tidal disruption 
events at the present epoch could be significantly more frequent than in models based on 
collisional loss-cone repopulation. The decay rate of a binary SBH could be signficantly 
enhanced compared with the rate in a spherical or axisymmetric nucleus, allowing binary 
SBHs to overcome the "final parsec problem" (Milosavljevic & Merritt 2003a) and coalesce 
by emission of gravitational radiation. Even if long-lived triaxiality should turn out to 
be rare, we show that transient departures from axisymmetry during mergers or galaxy 
interactions could induce feeding at rates approaching those inferred in AGNs. 



2. Models and Units 

The models of Papers II and III were based on the density law 

p* = p m~ 7 , (la) 
2 x 2 y 2 z 2 
a 2 b 2 c 2 

with 7 = (1, 2) and c/a = 0.5. The outer boundary (an equipotential surface) was chosen to 
contain roughly 100(20) times the black hole mass for 7 = 1(2). Non-evolving solutions for 
both values of 7 were found for T = 0.25 (nearly oblate) and T = 0.50 (maximal triaxiality); 
T = (a 2 — b 2 )/(a 2 — c 2 ) is the triaxiality index. Models with T = 0.75 (nearly prolate) evolved 
rapidly into precisely axisymmetric shapes and are not considered here. The black hole was 
represented by a central point with unit mass; the long-axis scale length a and the constant 
of gravitation G were also set to unity. The scale-free nature of the mass distribution allowed 
us to set p — 1 without loss of generality. 

Our focus is on models which include chaotic orbits. Four such models were presented 
in Paper III. Table 1 gives the axis ratios and the mass fractions on chaotic, tube and 
pyramid (box) orbits for these four models. (For definitions of the orbit families, see Papers 
I-III.) Chaotic orbits and tube orbits contributed roughly equally to the total mass in the 
self-consistent solutions, with regular box orbits a distant third. We therefore do not treat 
the box orbits separately from the chaotic orbits in what follows. The mass fraction in 
"centrophilic" orbits (chaotic orbits and pyramids) was ~ 50% in all models. 
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It is convenient to define equivalent spherical models for the stellar mass distribution; as 
shown in Paper III, these spherical models have energy distributions N(E) similar to those 
of the triaxial models and are useful when scaling the results found here to real galaxies. The 
equivalent spherical models have p*(r) = (r/5)~ 7 with 5 = (a&c) 1 / 3 , or 5 = 0.734(0.767), T = 
0.5(0.25). Their gravitational potentials in model units are 

$*(r) = 2n5r, 7 = 1 (2a) 



4ir5 2 



In [ ? ] - 1 



7 = 2. (2b) 



,5, 

The additive constants in these expressions for the potential have been chosen in the same 
way as for the triaxial models (see Paper I); we note that, for 7 = 2, the zero point of the 
potential occurs at a radius containing 4neS 3 times the black hole mass M., or 13. 51(15. 41)M. 
for T = 0.5(0.25). Note that the black hole has been omitted in these expressions for the 
potential. 

We define r\ to be the radius in the spherical model containing a mass in stars equal to 
twice the black hole mass (cf. Merritt 2003). This definition is equivalent to the standard 
one, r h = GM,/a 2 , when 7 = 2, the singular isothermal sphere. For 7 = 1, r h = (7rc))~ 1//2 = 
0.659(0.643) for T = 0.5(0.25), while for 7 = 2, r h = {2-K5 2 )- 1 = 0.296(0.270). The ID 
stellar velocity dispersion a is equal to \f2~nb = 1.84(1.92) in model units for 7 = 2. The 
velocity dispersion in a nucleus with 7 = 1 depends on radius and is zero at the center in 
the absence of a black hole (Dehnen 1993); hence we do not quote a value for a when 7 = 1. 

We define E h = $(r h ). In model units, E h = 3.04(3.10) for 7 = 1 and -12.93(-15.01) 
for 7 = 2. 

The results presented below can be related to real galaxies using the following scale 
factors for mass, length and time, derived from the spherical models just defined. Numerical 
constants, when given, are for T = 0.5 and differ only slightly for T = 0.25. 

[M] = M. (3) 



[L] = (n5f\ h ^ 152 pc , 7 = 1 (4a) 

= (2n5 2 ) r h * 37.4 pc ( T ) Y— 7 = 2 (4b) 

V ' y V200 kms-V V1O 8 .M ; [ V ; 



\n (^N/T^r-^xlOVf^V'V^V 172 , 7 = 1 (5a) 



gm. ' J Vioopcy V. 108 -^© 
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= (2^ 2 ) 3/2 \l4r = 3-30 x 10 5 yr ( ( J£-) , 7 = 2. (5b) 

V ; \ GM. y V200 kms-V V 10 Md/ 

In equation (5b), the relation a 2 = GM,/rh was used. 

The tidal disruption radius r t of a 10 8 A4 Q black hole is roughly equal to its Schwarzschild 
radius r s = 2GM,/c 2 . In model units, the Schwarzschild radius is 

— w 6.31 x 10 — 3— — , 7 = 1 (6a) 

1 CM. /cr\2 _ / a \ 2 

= ,o o ~ °- 591 - ~ S- 63 x 10~ 7 T , 7 = 2. 6b 

vr5 2 c 2 r h \cJ \200 km s^V V ; 

In order to keep the notation as simple as possible, the same symbols will be used for 
both dimensional and dimensionless quantities in what follows. When not noted explicitly, 
the distinction will be clear from the form of the equations. 

3. Pericenter Distributions 

Chaotic orbits pass near the central singularity once per crossing time, although the peri- 
center distance - the distance of closest approach to the black hole - varies quasi-randomly 
from passage to passage. We computed the statistics of pericenter passages by carrying out 
integrations of chaotic orbits for 10 5 To in each of the four triaxial potentials; T D is the 
energy-dependent dynamical time defined in Paper I. To speed the integrations, the poten- 
tial and forces due to the stars (equation lb) were expressed in terms of a truncated basis 
set (Hernquist & Ostriker 1992); Ten radial and angular functions were used. The routine 
RADAU (Hairer & Wanner 1996) was used for the integrations. Pericenter distances were 
computed by locating two times (ti, t-z) between which the distance to the black hole reached 
a minimum and interpolating the solution on a fine grid in ti < t < t 2 to find the distance 
of closest approach. 

Figure 1 shows the cumulative distribution of pericenter distances for two chaotic orbits 
integrated in the triaxial potentials with 7 = (1,2) and T = 0.5. The number of central 
encounters per unit time with pericenter distances less than d, N(r p < d), is nearly lin- 
ear with d over its entire range. The linear dependence was observed to extend down to 
pericenter distances of 10~ 6 or smaller, of order the tidal disruption radius in model units. 
The approximately linear dependence of N E on d, combined with the gravitational focussing 
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equation 

with Tb the impact parameter, implies N(r p < r&) oc r^, i. e. a nearly uniform filling of the 
two-dimensional cross section that defines the "throat" of the chaotic orbit (cf. Gerhard & 
Binney 1985). 

Chaotic orbits were observed to quickly fill the configuration-space region accessible to 
them, and one therefore expects that the statistical properties of one chaotic orbit at a given 
energy will be the same as those of any other chaotic orbit at the same energy. We verified for 
a number of chaotic orbits that N(r p < d) is indeed only weakly dependent on the starting 
point for fixed energy E. Hence we can write N(r p < d) = N E (r p < d), the rate of pericenter 
passages for a star on a chaotic orbit of energy E. 

We computed Ne{t p < d) at the energies associated with each of the mass shells defined 
in Paper II. Let A(E)d be the rate at which a single star on a chaotic orbit of energy E 
experiences pericenter passages with r p < d. The dependence of A on E in the four triaxial 
potentials is shown in Figure 2. For 7 = 2, an exponential fits the data well: 

InA^a + bE, (8) 

with 



a = -0.603, b = -0.290 (T = 0.5) (9a) 

a = -0.655, b = -0.275 (T = 0.25). (9b) 

We note that b —1/2tt5 2 . Henceforth we set b = —1/2tt5 2 ; thus A oc e~ E l° 2 and A scales 
with radius as ~ r~ 2 . This simple scaling will be useful in what follows. In physical units, 

A(S)«1.2^e-( M ^. (10) 

Figure 2 shows that this relation holds even at energies E < E h . 

For 7 = 1, a power law provides a good fit in the energy range E > E h : 

\nA^ a + b\nE, (11) 

with 



a = 1.734, b = -1.388 (T = 0.5) 
a = 1.869, b =-1.357 (T = 0.25). 



(12a) 
(12b) 
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In physical units, 



Thus A ~ r- 1A , r > r h . 



The pericenter distribution at any E is also characterized by a second quantity, the 
maximum pericenter distance r P)max [E) reached by chaotic orbits of energy E. This distance 
corresponds roughly to the width of the "throat" of a box orbit in an integrable triaxial 
potential. We find, for 7 = 2, 

\nr p , max pa c + dE, (14) 



with 



c = -0.43, d = 0.146 (T = 0.5) (15a) 
c = -0.36, d = 0.141 (T = 0.25). (15b) 



We note that d = 1/4tt5 2 to within its uncertainties; thus 



1 p,max ~ 



0.3r h e^ E ~ E ^ 2 (16) 



and r Pjmax scales approximately linearly with apocenter distance, r PtTnax pa 0.2r apo . 
For 7 = 1, we find, for E > E h , 

\nr Pimax ^ c + dlnE, (17) 

with 

c = -1.86, rf=0.86 (T = 0.5) (18a) 

c=-2.03, rf = 0.90 (T = 0.25). (18b) 



This implies 



g \ 0.88 



r p ,max ~ 0.60r A ( — ) (19) 



and again r p ^ max scales roughly linearly with r apo . 



4. Feeding Rates 



The feeding rates implied by Figure 2 are high compared with those in diffusive loss 
cone models; in fact they are comparable to the maximum possible capture rates in spherical 
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models, the so-called "full loss cone" rate. We illustrate this by considering the rate at which 
stars on a single orbit pass within a distance r t of the center. In the triaxial models, this 
rate is ~ A(E)r t . Setting E to the energy of a circular orbit at r/j, taking r t m 2 x 10~ 7 in 
model units (equation 6b for 7 = 2), and using the expressions given above for A(E), we 
find a rate of ~2x W^Tp 1 . In other words, roughly 10 6 orbital periods are required for 
a star on the edge of the black hole's sphere of influence to pass within its tidal radius. In 
a spherical model with a full loss cone, a fraction ~ r t r h /r 2 of the stars at r come within 
r t each orbital period; hence the mean capture rate from stars near is ~ (rt/r^T^ 1 or 
~ 1.5 x lCr 6 !^ 1 . In the spherical geometry, the net rate of pericenter passages is due to 
a small population of stars (those with sufficiently low angular momenta) being lost over a 
short period of time (an orbital period). In the triaxial geometry, a large population of stars 
(all stars on chaotic orbits) are supplied to the center over a longer period of time (~ 10 5 
periods). Of course, in the absence of loss cone repopulation, capture rates in the spherical 
geometry would drop to zero after a single orbital period. 

Next we compute the full feeding rates in the triaxial models. Let Ai c (E,t)dE be the 
mass in stars on chaotic orbits with energies from E to E + dE; the time dependence reflects 
a possible loss of stars due to tidal disruption or capture. If the capture or disruption radius 
is r t , the energy- dependent loss rate is 

M c (E,t)dE = —r t A(E)M c (E, t)dE (20a) 
= -r t A(E)M c (E, 0)e~ A(E)nt dE (20b) 

and the total capture rate from chaotic orbits at all energies is 

M(t) =r t J A(E)M C (E, 0)e~ A{E)nt dE. (21) 

We evaluated these expressions in two ways. Ai c (E, 0) can be computed directly from 
the orbital weights in the Schwarzschild solutions (cf. Figure 3 of Paper III). Alternatively, 
smooth approximations to A4 C (E, 0) can be computed using the equivalent spherical models 
defined above. In the latter case, we define M C {E, 0) = f c (E)M{E) with M{E) the energy 
distribution in the equivalent spherical model and f c (E) the fraction of stars at energy E 
assumed to be on chaotic orbits. For 7 = 1, the energy distribution in the spherical geometry 
is 



or, in model units, 



M( V = T 5 gW. E (22a) 



M(E) = —E. (22b) 
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For 7 = 2, we have 



M(E) 



2y/Er h 
9 G 



e (E-E h )/2a 2 



(23a) 



or, in model units, 



M(E) 



2V6 



eSe E ^ s \ 



(23b) 



9 



Since these approximations to A4(E) were derived from expressions for the potential that 
exclude the central point mass, they are only strictly valid outside the black hole's sphere 
of influence, E > Eh- However the expression for 7 = 2 turns out to be reasonably correct 
even at lower energies (cf. Figure 3 of Paper III), a result that will be used below. 



set to a constant, f c , in each model; no attempt was made to match the detailed energy 
dependence of the chaotic orbital populations in the Schwarzschild solutions. Nevertheless 
the fits are reasonably good, particularly for 7 = 2. For 7 = 1, the orbital distributions in 
the Schwarzschild solutions are "noisy" but the analytic expression does a reasonable job 
of reproducing the mean dependence. The values of f c used in Figure 3 were f c = 0.25 
( 7 = 2,T = 0.5 and 0.25), J c = 0.5 (7 = 1,T = 0.5) and J c = 0.4 (7 = 1,T = 0.25). These 
values are slightly smaller than the overall chaotic mass fractions in the numerical solutions 
(Table 1); the reason is that the Schwarzschild solutions put a large number of orbits in the 
outermost shells to compensate for the hard outer edge, and are correspondingly depleted 
at low energies (Figure 3 of Paper III). 

Also shown in Figure 3 are the feeding rates predicted by a spherical, "full loss cone" 
model. The full-loss-cone feeding rate is 



where M.i c (E) is the number of stars in the spherical model with pericenters below r t , 
f(E) is the isotropic distribution function, P(E) is the radial period, and J/ c is the angular 
momentum of a star with pericenter 77, Jf c = 2r\{E — $(r t )) ps 2GM,r t . Equation (24b) 
describes the capture rate in a spherical model assuming that all of the orbits remain fully 
populated. Figure 3 shows that the capture rates in the triaxial models are similar to the 
full-loss-cone rates in the equivalent spherical models, and can even exceed them at high 
energies. 



Figure 3 shows Ai c (E, 0) computed in both ways. The chaotic orbit fraction f c (E) was 



M full (E)dE 




(24a) 



A^f{E)Jl{E)dE 



(24b) 



In a real spherical galaxy, the high feeding rates corresponding to a full loss cone would 
persist for just a single orbital period at each energy. In the triaxial models, by contrast, 
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the capture rate from chaotic orbits decreases only by a factor ~ exp [A(E)P(E)r t ] in one 
orbital period. At energies > E h , this factor is negligible meaning that feeding rates will 
remain large for many orbital periods. 

The evolution of M c (E,t) is shown in Figure 4. Here we have set the capture radius 
r t to its smallest possible value, r t = r s , using equation (6b). For 7 = 1, this requires a 
choice for M. and r^; we took M. = 1O 8 A^0 and = 100 pc. For 7 = 2, only a needs to 
be specified; we took a = 200 km s _1 . The analytic approximations, equations (22b) and 
(23b), were used for M c (E,0) and we set f c = 0.5. Figure 4 shows that, for 7 = 1, the 
capture rate is low enough that no significant changes occur in Ai c (E, t) even at times as 
great as t — 10 6 , corresponding to ~ 10 12 yr using the scaling of equation (5b). For 7 = 2, 
the capture rate is higher and changes in Ai c {E, t) begin to occur at energies above E h for 
t > 10 5 , corresponding to ~ 10 10 yr. 

The total capture rate from orbits at all energies is M(t) = f M. c (E,t)dE. For 7 = 1 
and r t r s , the captured mass is a small fraction of M. and rh is nearly independent of 
time. Thus we can write 

-Ei r E 2 



M « / M c (E,0)dE = r t A(E)f c (E)M(E)dE (25a) 

J E\ J E\ 

rE2 

z a r t / f c (E)E 1+b dE (25b) 

J Ei 



35tt8 j El 

in model units. As lower integration limit we take E\ = 0; this is reasonable since the value 
of the integral is not strongly dependent on the lower cutoff. However since 1 + 6 ~ —0.4, the 
integral diverges for large E 2 . This corresponds physically to the fact that a p oc r _1 nucleus 
can not extend indefinitely. Nuclear density profiles in real, weak-cusp galaxies exhibit a 
break at a radius r b IOtv We define E b = $(r b ) and set E 2 = E b) giving, in model units, 

M = , 8 - e a r t J c E 2 h +h . (26) 
35(2 + 6)tt5 tJc b K J 

Setting r t — r s , this becomes, in physical units, 

M =(1.77, 1.97) xlO-^ yr-^ -^- (27) 



100 pc/ \10 S M Q J \10r h 

where the number pairs in parentheses refer to T = (0.5, 0.25) respectively. We note that 
the feeding rate is almost the same for the two values of the triaxiality index, suggesting a 
weak dependence of M. on the degree of departure from axisymmetry, for a given chaotic 
mass fraction. The feeding rates of equation (27) imply a captured mass of < 0.01M. over 
a Hubble time. 
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For 7 = 2, the feeding rates are much higher, and the captured mass can be a large 
fraction of M. even when r t is as small as r s . Hence the population of chaotic orbits is 
significantly depleted and we can not ignore the time dependence of m. We have 

M(t) = r t J A(E)M c (E,0)e- A(E)rtt dE (28a) 

= ^^e a r t Jf c (E)e- E ^ s2 exp(r t te a+bE )dE (28b) 

in model units. It is reasonable to set the lower integration limit to — oo, since our expressions 
for A(E) and M.(E, 0) are valid for E <C E h and the contribution to the integral from 
low energies falls off quickly with time as the most-bound particles are eaten; hence any 
errors due to the forms of A(E) or Ai(E) at low energies are transient. Setting the upper 
integration limit to +oo is also reasonable since the product A(E)At(E, 0) drops rapidly 
with E. (Furthermore it is possible that the combined, luminous plus dark matter density 
profiles in early- type galaxies are well described as p ~ r~ 2 even far outside of the nucleus.) 
Making the substitution y = r t texp(a + bE), we find 

M (t) = W6^>w^ i/2 r r ^ e -» dy (29a) 



1 /2 

» (4.82, 5.36) 7c (j) (29b) 

where the paired numbers refer to T = (0.5,0.25). We note again the small dependence of 
the feeding rate on the degree of triaxiality for fixed f c . The total mass accreted after time 

tacc IS 

AM = (9.64, 10.72) J c {r t t acc ) l/2 . (30) 
In physical units, these expressions become 

3 / \ 1/2 / . \ -1/2 



M(t) = (1.48, 1.54) { (31a) 



(4.28,4.46) x lO- 3 ^ yr-V c ^ 

/ M. \ 1/2 / t N - 1/2 
\10 8 M Q ) \10 10 jt 



1/2 / ^ \ 5/2 

200 km s _1 ' 



(31b) 



and 



/ \ 1/2 , . \ 1/2 

AM = (2.97, 3.09) TM."-^) (^) (32a) 
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.56,8.92) x!O.M /^-j ^^rj x 
f^Y'Vj^Y' 2 (32b) 



V.io 8 .M y V!0 10 yr 

1/2 

The r t dependence reflects the more rapid destruction of the cusp when r t is large. Equa- 
tion (32b) implies a captured mass of order M. in a Hubble time. 



5. Tidal Disruption Rates 

The feeding rate in a steep-cusp nucleus is given by equation (31b), as a function of the 
capture radius r t . If M. < 1O 8 A^ , solar type stars are tidally disrupted at a radius 

/ M. V 2/3 

(Hills 1975). (10 8 .M© is also roughly the mass of the largest black holes that sit in steep- 
cusp nuclei.) Combining equations (31b) and (33), the tidal disruption rate in a steep-cusp 
nucleus is 

M colllsimless (t) « 4 x 10' 3 M Q yr^Tc a^M^t'^ 2 (34) 

with tio the time since cusp formation in units of 10 10 yr, <t 2 oo the velocity dispersion in 
units of 200 km s" 1 and M. j8 the black hole mass in units of 10 S M Q . We refer to this as the 
"collisionless" tidal disruption rate to highlight that the supply of stars to the black hole is 
not being driven by gravitational scattering, as in the standard model (Frank & Rees 1976; 
Lightman & Shapiro 1977). 

Supermassive black holes in galaxies at the current epoch satisfy the M. — a relation, 



f M. \ ....,/ a 



4.65 



, . . .« 1.48 t (35) 

VlO 8 A^ y V200 kms~V V ; 

(e.g. Merritt & Ferrarese 2001b), which allows us to write the tidal disruption rate in nearby 
galaxies in terms of either M. or a alone: 

M colllsionless (t) « 5 x 10-^Qyr- 1 ^ t^ /2 (36a) 
» 4 x lO^eyr- 1 ^ M°J° t^ 2 . (36b) 
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We compare this to the rate at which collisional (encounter-driven) loss-cone refilling 
supplies stars to the central black hole. In a spherical, p oc r~ 2 nucleus, 

M coUlslonal = 2.8 x KT^oyr- 1 <t 7 £m^ (37) 

(Wang & Merritt 2003); this expression again assumes that the disrupted stars have the mass 
and radius of the sun. Equation (37) can likewise be rewritten using the M. — a relation: 

Mcoiusionai ~ 2 x 10' 4 M yr~ V^- 15 (38a) 
w 2 x lO^AlQyr^M-g - 25 . (38b) 

Figure 5 plots M(M.) in the collisionless (triaxial) and collisional (spherical) cases for a 
p ~ r~ 2 nucleus. The two rates scale in the opposite sense with M. and collisional feeding 
would dominate in sufficiently small galaxies. The rates are equal when 

M. w 5 x lO^oTT 1 ^ 5 - (39) 

1 /2 — 

For t 10 / f c pa 1, this mass is close to that of the smallest black holes with reliably determined 
masses, in the Milky Way and M32. Between this mass and M. pa 1O 8 A^0, Figure 5 suggests 
that loss cone feeding driven by triaxiality can easily dominate collisional loss cone feeding, 
as long as the chaotic mass fraction f c is not too much smaller than one. Flaring rates should 
peak in the brightest galaxies with steep nuclear density profiles, at values of a few times 
10~ 3 yr _1 . Since the disruption rate in the triaxial geometry varies as t~ 5 oc (1 + z) - 5 , 
flaring could be even more important at intermediate redshifts. 

In a p oc r _1 nucleus, and again assuming M, < 1O 8 A^ , equation (27) gives a disruption 
rate for solar-type stars of 

_ / r , \ - 5 / 2 / M \ 11/6 / r h \ 5/8 
M coUlslonless ^2xlO- 5 M Qj T- 1 fJ—^) -T^H 777^ > ( 40 ) 



100 pcy \10 S M Q J \10r h 
approximately independent of time. The corresponding collisional expression is 

_ / r , \ -V3 / M \ 29/18 
M collisional » 1 x lO^ 6 ^ yr-V c — h — . (41) 



vioopcy \io 8 m q/ 

The functional dependence in equation (41) is taken from Wang & Merritt (2003), and the 
normalizing factor is based on those authors' calculation of iV for the galaxy NGC 3379 
(7 pa 1.1). For a 7 = 1 nucleus, the collisional and collisionless disruption rates scale in 
almost the same way with M. and r^, and the collisionless rate exceeds the collisional rate 
as long as f c > 0.05. 
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In summary: even small fractional populations of chaotic orbits can produce tidal dis- 
ruption rates that exceed those predicted by the standard model of collisional loss cone 
repopulation in a spherical nucleus (Syer & Ulmer 1999; Magorrian & Tremaine 1999; Wang 
& Merritt 2003). In triaxial nuclei with p ~ r -2 and M. m 1O 8 A^0, tidal flaring rates can 
plausibly reach values as high as several times 10~ 3 yr _1 for solar type stars. Rates could be 
even higher if the cusps were recently formed. 



6. Black Hole Growth and the M, — a Relation 

The fate of gas liberated by the tidal disruption of a star at r t > r s is uncertain. But 
for M, > 1O 8 .M , stars are directly captured by the black hole without being disrupted. 
Using equation (32b), setting r t = r s and requiring AM = M., we find the following relation 
between M. and a in a steep triaxial cusp: 

M - mJJ^f^JL^X. ,4 2 ) 



10 S M Q ' JC \10 10 yr J V200 km s" 1 

This is remarkably similar to the M. — a relation: the exponent on a is consistent with 
measured values, 4.5 ± 0.5 (Merritt & Ferrarese 2001a) and even the normalization is of the 
right order if f c w 1. 

The M oc <t 5 dependence can be understood in several ways. Writing M.(E, 0) = 
r s A(E)M c (E ) 0) and integrating over energies from E h to infinity, we find 

M(E> E h )^ 2.5 /c^(£) 2 • (43) 

Alternatively, the total rate at which stars pass inward at radius r is ~ A-Kr 2 p(r)a = 2<r 3 /G, 
and a fraction F = Tr>Af c r s of stars at energy E pass within r s each crossing time. Equation 
(10) gives 

a r 



Air) *GHl?{7 h ) (44) 
and T D oc (r h /a)(r/r h ), so that 

M(r)0C ^© 2 7- (45) 

These relations show that the feeding rate due to stars originating outside of the SBH's 
sphere of influence is of order 
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independent of M. - large enough to contribute substantially to the masses of SBHs if 
accretion continues for 10 9 yr or more. 

For M. < 10 8 M & , the tidal disruption radius exceeds r s and scales as M. 1 ^ 3 (equation 
33). If the growth of the black hole is determined by the rate of tidal disruptions, AM scales 
more weakly with a in this regime, roughly as a 3 . This would imply a flattening of the M. — cr 
relation at masses below ~ 10 8 .A/f©. However it is not clear what the fate of tidally-liberated 
gas would be; some would fall into the black hole but some would escape (Frank 1979). 
For masses below ~ 10 6 A4 Q , Figure 5 suggests that collisional loss cone repopulation would 
dominate over collisionless feeding. 

Equation (46) is similar to equation (8) of Zhao, Haehnelt & Rees (2002), who considered 
the consequences for black hole growth of a continuously resupplied loss cone in the spherical 
geometry. The agreement is reasonable since the feeding rates in the triaxial models are 
comparable to spherical, full- loss-cone rates for E > Eh (Figure 3). The triaxial models are 
different in one respect. In the spherical geometry, the contribution per orbital period to 
AM from stars at radius r scales as r~ 2 . In the triaxial models, every centrophilic orbit is 
capable of visiting the center and the probability that a single star will do so in one orbital 
period varies as ~ r -1 . Thus while the feeding rates from stars at r r\ are comparable 
in the two models, large-radius orbits contribute relatively more in the triaxial models. One 
consequence is that the growth of the black hole is not so strongly dependent on the stellar 
distribution near r h . 

Accretion of stars, or of gas liberated from stars, have long been discussed as mechanisms 
for growing SBHs in galactic nuclei (Hills 1975; Young, Shields & Wheeler 1977; Frank 1978; 
McMillan, Lightman & Cohn 1981). Hills pointed out already in 1975 that feeding at the 
full-loss-cone rate could plausibly grow a 1O 8 A^ black hole in a Hubble time. Subsequent 
authors (e.g. Young, Shields & Wheeler) tended criticize this claim on the grounds that 
gravitational encounters occur too infrequently to maintain a fully-populated loss cone. Our 
result - that sustained feeding rates in the triaxial geometry can be comparable to full-loss- 
cone rates in the spherical geometry - should help to revive Hills' model. However, chaotic 
loss cones still fail to reproduce in a natural way the observed time dependence of quasar 
luminosities. The t~ 1//2 time dependence found here for stellar feeding in a p ~ r~ 2 density 
cusp (equation 31b) - which results entirely from the depopulation of chaotic orbits - is 
much more gradual than the observed, ~ (1 + z) 3 falloff of quasar luminosities (e.g. Boyle, 
Shanks & Peterson 1988), and would require that much of the growth took place long after 
the quasar epoch. Slow, optically faint growth of SBHs has been proposed (Richstone et 
al. 1998), but only before it was discovered that SBH masses in nearby galaxies had been 
substantially overestimated (Merritt & Ferrarese 2001b). SBH masses are now consistent 
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with accretion rates inferred from quasar statistics (Merritt & Ferrarese 2001c) leaving less 
room for growth at late times. 

On the other hand, the time dependence derived here assumes that the steep density 
cusp forms only once and that it is slowly depleted by accretion of centrophilic orbits. Cusps 
might form episodically via star formation during mergers leading to repeated bursts of 
stellar fueling and larger average rates of accretion (Milosavljevic & Merritt 2003b). Other 
considerations could boost the mean capture rate as well. Growth of the SBH via accretion 
steepens the stellar density profile, to p ~ r ~ 2 - 5 near the SBH (Merritt 2003), an effect that 
was ignored here. Depopulation of the chaotic orbits would be counteracted to some extent 
by scattering of stars onto these orbits. Dark matter could also be accreted (Zhao, Haehnelt 
& Rees 2002; Read & Gilmore 2003). Once the supply of star-forming gas was depleted 
(Kauffmann & Haehnelt 2000), subsequent mergers would turn off the stellar feeding through 
the formation of a binary SBH which ejects stars from the nucleus and lowers its density 
(Milosavljevic & Merritt 2001). This might be a natural way to explain the rapid observed 
falloff in quasar luminosities. In summary, while capture of stars is probably not the only way 
in which black holes grow, there would seem to be no reason to rule it out as an important 
contributor. 



7. Transient Changes in the Potential 

Deviations from axial symmetry are likely to occur at least temporarily following a 
merger or accretion event, yielding chaotic orbits near the black hole and resulting in en- 
hanced feeding rates (Norman & Silk 1983). Assuming a steep cusp (perhaps formed dissi- 
patively during the merger), equation (32b) and the M. — a relation imply an accreted mass 
of 

/ \ 0.5 / \ 5 / A , \ 0.5 

ite over At 

At oy \r s J V200 kms-V \10 8 yr/ 

This is comparable to feeding rates inferred in AGNs for a < 250 km s _1 and in quasars for 
a > 350 km s" 1 . 



AM w 2 x 10 7 M Q f c - T — — (47) 

1 ' \2ookms-V v 108 y r , 



in time At, or a mean feeding rate over At of 
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8. Decay of a Black-Hole Binary 

Following the merger of two galaxies each containing a SBH, a binary forms with semi- 
major axis a = dhard, where 

^ * ^ W 2 ' 7 PC(1 + ^ (A) fc) 2 (49a) 
" '•/,: (49b) 



4(1 +p) 

m2 is the mass of the smaller black hole, p = 1712/1711 and /i = m^/M, (Merritt 2003). 
In a spherical or axisymmetric galaxy, the binary quickly ejects stars with pericenters below 
~ cihard, after which changes in the binary separation take place on the much longer time 
scale associated with collisional refilling of the binary's loss cone (Yu 2002). In all but the 
densest nuclei, decay stalls at a separation a > 0.1a hard , too large for the efficient emission 
of gravitational waves (Merritt 2003). In a triaxial galaxy, decay can continue due to the 
much larger supply of centrophilic stars. The decay rate is determined by the supply of stars 
into a region of radius ~ a around the center of mass of the binary (Quinlan 1996). In a 
steep triaxial cusp, the feeding rate from stars of energy E > &(ah ar d) into a region of radius 
~ dhard is given by an equation similar to (43): 

POO 

M(E > E hard ) « a hard / A(E)M c (E)dE (50a) 

J Ehard 

» 1.3/^ (50b) 

« 25OOM yr- 1 /; ( °- ^ , (50c) 

V200 kms 3 y V ; 

independent of M. and p. If this rate were maintained, the binary would interact with its 
own mass in stars in a time of only ~ 10 5 yr. In fact, the feeding rate will decline with time 
as the centrophilic orbits are depleted. Equating the energy carried away by stars with the 
change in the binary's binding energy gives 

-dM ss ^-^d - (51) 



2 a 2 \a, 

(Merritt 2003). The coupled equations describing the change in the binary separation and 
the evolution of the stellar distribution are then 



d(l/a) 
dt 

dM c 
dt 



= All f A(E)M c (E,t)dE, (52a) 
M. a J 

= —r t A(E)M c (E,t) (52b) 
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where r t is the pericenter distance below which stars can efficiently exchange energy with 
the binary. Note that we are ignoring any contribution to shrinkage of the binary from 
low-angular-momentum regular orbits, or from collisional repopulation of the loss cone. 

We set r t to a (small) multiple of a, r t (t) = 1Za(t), 1Z « 1, reflecting the fact that 
only stars which come to within roughly a binary separation of the black holes will inter- 
act strongly enough with them to be ejected. As lower integration limit on the integral in 
equation (52a), we take E^ard] this choice reflects the fact that the stellar distribution imme- 
diately following the formation of a hard binary is similar to that in the pre-merger galaxies 
outside of r pa ahard (Milosavljevic & Merritt 2001). The evolution equations become 

d(l/a) 3 

' Ehard 

= -Ka(t)A(E)M c (E,t). (53b) 



dt 

dM c 



3 f°° 

—11 / A(E)M c (E,t)dE, (53a) 

M » JE hn „, 



dt 



As initial conditions we take a(0) = ahard and Ai c (E, 0) = f c A4(E) with A4(E) given by 
equation (23b). 

Figure 6 shows solutions for various choices of f c and for p = m 2 /mi = (1, 0.1). At late 
times, the binary separation varies approximately as 

^nffxt (54) 

in model units, for both values of p. In physical units, 

^f- « 3 x 10 4 ^ 2 (2oo]^0 (l0^) (lO^v)' (55) 

A 1/a oc t dependence is expected in "full loss cone" situations like this one (Milosavljevic 
& Merritt 2003b). The nonlinear dependence of a(t) on f c reflects the fact that, when f c is 
small, the decay rate is also small and the binary tends to stall at a large separation, causing 
it to interact with and deplete the supply of stars, further slowing the decay etc. 

Coalescence due to gravitational wave emission in a time t gr occurs when a = a gr , where 

^ « 75^ ( * \ ~ 7/8 (Jl-) ~ 1/4 (56) 
a gr (1+pf 2 V200 kms-V V10 9 yr/ 

for a circular-orbit binary (Merritt 2003). For p = 1(0.1), equation (56) implies that decay 
by a factor of ~ 10 2 (10 1 ) is required in order for gravity-wave coalescence to occur in 10 9 
yr. Combining equations (55) and (56), this can be achieved in 10 10 yr for chaotic mass 
fractions: 

_ / a \9/i6 

005 (siiotoFr) (57) 
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with a weak dependence on p. Thus, placing just a few percent of a galaxy's mass on chaotic 
orbits is sufficient to overcome the "final parsec problem" (Milosavljevic & Merritt 2003a) 
and induce coalesence, if (as assumed in Figure 6) the chaotic orbits are present initially at 
energies much greater than E^. Since the stars supplied to the binary come from well outside 
of Th in this case, the effect on the stellar density profile would be negligible even after the 
binary had ejected several times its own mass in stars. 



9. Summary 

In a p ~ r~ 2 triaxial nucleus containing chaotic (centrophilic) orbits and a central black 
hole, stars are supplied to the black hole's sphere of influence at a rate ~ f c a 3 /G, and to the 
black hole itself at a rate ~ f c cr 5 /Gc 2 , where f c is the fraction of the mass on centrophilic 
orbits and a is the stellar velocity dispersion. If f c is of order unity, as in recently published 
self-consistent models (Papers II, III), the mass accumulated in 10 10 yr is of the same order 
as observed black hole masses. Feeding rates in such an environment fall off as ~ t^ 1 / 2 as 
the chaotic orbits are depleted. Tidal disruption rates at the current epoch due to stars on 
chaotic orbits are as much as 1-2 orders of magnitude greater than in standard, spherical 
or axisymmetric models in which loss-cone refilling occurs via gravitational scattering onto 
eccentric orbits. Transient changes in the shape of a nucleus during mergers or accretion 
events could lead to episodic stellar feeding at rates of ~ 0.2 — 2Ai Q yr -1 , comparable 
to the rates inferred in active galactic nuclei and quasars. Decay of a black-hole binary 
at the center of a triaxial nucleus could be greatly enhanced compared with spherical or 
axisymmetric nuclei, even if only a few percent of the mass is on chaotic orbits, thus solving 
the "final parsec problem." 

This work was supported by NSF grants AST 00-71099 and AST 02-0631, and by NASA 
grants NAG5-6037 and NAG5-9046. M.Y. Poon is grateful to the Croucher Foundation for 
a postdoctoral fellowship. 



-20- 



REFERENCES 

Boyle, B. J., Shanks, T. & Peterson, B. A. 1988, MNRAS, 235, 935 

Crane, P. et al. 1993, AJ, 106, 1371 

Dehnen, W. 1993, MNRAS, 265, 250 

Dubinski, J. & Carlberg, R. 1991, ApJ, 378, 496 

Erwin, P. & Sparke, L. S. 2002, AJ, 124, 65 

Frank, J. 1978, MNRAS, 184, 87 

Frank, J. 1979, MNRAS, 187, 883 

Frank, J. & Rees, M. J. 1976, MNRAS, 178, 633 

Gebhardt, K. et al. 1996, AJ, 112, 105 

Gerhard, O. E. & Binney, J. 1985, MNRAS, 216, 467 

Hairer, E., & Wanner, G. 1996, Solving Ordinary Differential Equations II. (Berlin : 
Springer) 

Hernquist, L. & Ostriker, J. 1992, ApJ, 386, 375 
Hills, J. G. 1975, Nature, 254, 295 

Hughes, S. A., Marka, Szabolcs, Bender, P. L. & Hogan, C. J. 2001, eConf C010630, 402 
(astro-ph/01 10349) 

Kauffmann, G. & Haehnelt, M. 2000, MNRAS, 311, 576 

Lightman, A. & Shapiro, S. 1977, ApJ, 211, 244 

Magorrian, J. & Tremaine, S. 1999, MNRAS, 309, 447 

McMillan, S. L. W., Lightman, A. P. & Cohn, H. 1981, ApJ, 251, 436 

Merritt, D. 2003, to appear in Carnegie Observatories Astrophysics Series, Vol. 1, "Coevolu- 
tion of Black Holes and Galaxies," ed. L. C. Ho (Cambridge: Cambridge Univ. Press) 
(astro-ph/0301257) 

Merritt, D. & Ferrarese, L. 2001a, ApJ, 547, 140 



- 21 - 



Merritt, D. & Ferrarese, L. 2001b, in ASP Conf. Ser. Vol. 249, The Central Kiloparsec of 
Starbursts and AGN, ed. H. Knapen, J. E. Beckman, I Shlosman, & T. J. Mahoney 
(San Francisco: ASP), 335 

Merritt, D. & Ferrarese, L. 2001c, MNRAS, 320, L30 

Milosavljevic, M. & Merritt, D. 2001, ApJ, 563, 34 

Milosavljevic, M. & Merritt, D. 2003a, to appear in Proceedings of the 4th LISA Symposium 
(astro-ph/0212270) 

Milosavljevic, M. & Merritt, D. 2003b, ApJ, submitted (astro-ph/0212459) 

Norman, C. & Silk, J. 1983, ApJ, 266, 502 

Peng, C. Y., Ho, L. C, Impey, C. D. & Rix, H.-W. 2002, AJ, 124, 266 

Poon, M.Y. & Merritt, D. 2001, ApJ, 549, 192 (Paper I) 

Poon, M.Y. & Merritt, D. 2002, ApJ, 568, L89 (Paper II) 

Poon, M. Y. & Merritt, D. 2003, ApJ, in press (astro-ph/0212581) (Paper III) 

Quinlan, G. D. 1996, NewA, 1, 35 

Read, J. I. & Gilmore, G. 2003, MNRAS, in press (astro-ph/0210658) 

Rest, A. et al. 2001, AJ, 121, 2431 

Richstone, D. O. et al. 1998, Nature, 395A, 14 

Shlosman, I., Begelman, M. C. & Frank, J. 1990, Nature, 345, 21 

Syer, D. & Ulmer, A. 1999, MNRAS, 306, 35 

Valluri, M. & Merritt, D. 1998, ApJ, 506, 686 

Wang, J. & Merritt, D. 2003, preprint (astro-ph/0305493) 

Wozniak, H., Friedli, D., Martinet, L., Martin, P. & Bratschi, P. 1995, AA Suppl, 111, 115 
Young, P. J., Shields, G. A. & Wheeler, J. C. 1977, ApJ, 212, 367 
Yu, Q. 2002, MNRAS, 331, 935 

Zhao, H.-S., Haehnelt, M. & Rees, M. J. 2002, New Astron., 7, 385 



-22 - 




Fig. 1. — Number of central encounters per unit time with pericenter distance less than d, 
for two chaotic orbits, each integrated for 10 5 dynamical times in their respective potentials. 
The solid lines have unit slope. 



This preprint was prepared with the AAS IATgX macros v5.0. 
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Fig. 2. — The function A(E) that describes the cumulative rate of pericenter passages in 
the four triaxial models. Points are from integrations of chaotic orbits; lines show the fits 
described in the text. Arrows indicate Eh. 
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Fig. 3. — Energy-dependent capture rate Ai c (E) due to chaotic orbits in the four triaxial 
models ait — 0. The capture radius r t has been set to unity; capture rates scale linearly with 
r t (cf. eq. 20b). Open circles: Ai c computed using the actual chaotic orbit populations in 
the Schwarzschild models. Solid lines: M. c computed using the analytic approximations to 
M. c described in the text. Dashed lines: Capture rates predicted by a spherical, full-loss-cone 
model. Arrows indicate E h . 
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Fig. A.— M c (E,t) at four times, t = (0, 10 3 , 10 4 , 10 5 , 10 6 ) in model units, for 7 = 1 and 2, 
and T = 0.5. The capture radius r t has been set to the Schwarzschild radius r s using the 
scalings discussed in the text. Arrows indicate Eh- 
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Fig. 5. — Tidal disruption rates for solar type stars in collisional (spherical) and collisionless 
(triaxial) nuclei with p oc r~ 2 . Line labelled "collisional" is equation (37). "Collisionless" 
lines are equation (35) with f c ti^ 2 = (0.3, 1, 3). 
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Fig. 6. — Decay of a massive binary black hole at the center of a steep-cusp (p ~ r~ 2 ) nucleus 
due to ejection of stars on chaotic orbits, (a) m 2 /mi = 1; (b) m 2 /mi = 0.1. The different 
curves are for chaotic mass fractions of f c — 1 (top), 0.3,0.1,0.03 and 0.01 (bottom). Time 
is measured in model units (equation 5b); t = 10 5 corresponds roughly to 10 10 yr. 



